DJIA,GSPC,NDX,GDAXI,FCHI,SSEC,SENSEX Analysis

https://github.com/cran/VineCopula - Copula family

setwd('~/Documents/VaR')

prices <- read.csv('./data/StockIndexData.csv')

library(VineCopula)
library(copula)
## 
## Attaching package: 'copula'
## The following object is masked from 'package:VineCopula':
## 
##     fitCopula
library(tseries)

DJIA,GSPC,NDX,GDAXI,FCHI,SSEC,SENSEX

prices.DJIA <- prices[,c('DJIA')]

summary(prices.DJIA)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6547   10220   11020   12020   13270   18640
n <- length(prices.DJIA)
rtn.DJIA <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.DJIA[i] <- log(prices.DJIA[i+1]/prices.DJIA[i]) * 100.
}

prices.GSPC <- prices[,c('GSPC')]

summary(prices.GSPC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   676.5  1133.0  1296.0  1364.0  1478.0  2190.0
n <- length(prices.GSPC)
rtn.GSPC <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.GSPC[i] <- log(prices.GSPC[i+1]/prices.GSPC[i]) * 100.
}

prices.NDX <- prices[,c('NDX')]

summary(prices.NDX)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   938.5  1569.0  1944.0  2394.0  2954.0  4910.0    1063
n <- length(prices.NDX)
rtn.NDX <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.NDX[i] <- log(prices.NDX[i+1]/prices.NDX[i]) * 100.
}

prices.GDAXI <- prices[,c('GDAXI')]

summary(prices.GDAXI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2203    4934    6172    6461    7583   12370      60
n <- length(prices.GDAXI)
rtn.GDAXI <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.GDAXI[i] <- log(prices.GDAXI[i+1]/prices.GDAXI[i]) * 100.
}

prices.FCHI <- prices[,c('FCHI')]

summary(prices.FCHI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2403    3649    4220    4291    4836    6857      50
n <- length(prices.FCHI)
rtn.FCHI <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.FCHI[i] <- log(prices.FCHI[i+1]/prices.FCHI[i]) * 100.
}

prices.SSEC <- prices[,c('SSEC')]

summary(prices.SSEC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1011    1572    2104    2288    2850    6092     317
n <- length(prices.SSEC)
rtn.SSEC <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.SSEC[i] <- log(prices.SSEC[i+1]/prices.SSEC[i]) * 100.
}

prices.SENSEX <- prices[,c('SENSEX')]

summary(prices.SENSEX)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2600    4757   13670   13040   18840   29680     202
n <- length(prices.SENSEX)
rtn.SENSEX <- rep(0,n-1)

for(i in seq(n-1)){
  rtn.SENSEX[i] <- log(prices.SENSEX[i+1]/prices.SENSEX[i]) * 100.
}
rtn.DJIA.ar <- ar(rtn.DJIA, order.max = 1, method = 'ols')
rtn.DJIA.ar.resid.garch <- garch(rtn.DJIA.ar$resid[2:length(rtn.DJIA.ar$resid)],trace=FALSE)

rtn.GSPC.ar <- ar(rtn.GSPC, order.max = 1, method = 'ols')
rtn.GSPC.ar.resid.garch <- garch(rtn.GSPC.ar$resid[2:length(rtn.GSPC.ar$resid)],trace=FALSE)

rtn.NDX.1<-na.omit(rtn.NDX)
rtn.NDX.ar <- ar(rtn.NDX.1, order.max = 1, method = 'ols')
rtn.NDX.ar.resid.garch <- garch(rtn.NDX.ar$resid[2:length(rtn.NDX.ar$resid)],trace=FALSE)

rtn.GDAXI.1<-na.omit(rtn.GDAXI)
rtn.GDAXI.ar <- ar(rtn.GDAXI.1, order.max = 1, method = 'ols')
rtn.GDAXI.ar.resid.garch <- garch(rtn.GDAXI.ar$resid[2:length(rtn.GDAXI.ar$resid)],trace=FALSE)

rtn.FCHI.1<-na.omit(rtn.FCHI)
rtn.FCHI.ar <- ar(rtn.FCHI.1, order.max = 1, method = 'ols')
rtn.FCHI.ar.resid.garch <- garch(rtn.FCHI.ar$resid[2:length(rtn.FCHI.ar$resid)],trace=FALSE)

rtn.SSEC.1<-na.omit(rtn.SSEC)
rtn.SSEC.ar <- ar(rtn.SSEC.1, order.max = 1, method = 'ols')
rtn.SSEC.ar.resid.garch <- garch(rtn.SSEC.ar$resid[2:length(rtn.SSEC.ar$resid)],trace=FALSE)

rtn.SENSEX.1<-na.omit(rtn.SENSEX)
rtn.SENSEX.ar <- ar(rtn.SENSEX.1, order.max = 1, method = 'ols')
rtn.SENSEX.ar.resid.garch <- garch(rtn.SENSEX.ar$resid[2:length(rtn.SENSEX.ar$resid)],trace=FALSE)

rm(prices.DJIA, prices.GSPC, prices.NDX, prices.GDAXI, prices.FCHI,prices.SSEC, prices.SENSEX)
rm(rtn.DJIA, rtn.GSPC, rtn.NDX, rtn.GDAXI, rtn.FCHI,rtn.SSEC, rtn.SENSEX)
rm(rtn.NDX.1, rtn.GDAXI.1, rtn.FCHI.1,rtn.SSEC.1, rtn.SENSEX.1)
std.residual.DJIA <- rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.DJIA.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GSPC <- rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GSPC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.NDX <- rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.NDX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GDAXI <- rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GDAXI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.FCHI <- rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.FCHI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SSEC <- rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SSEC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SENSEX <- rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SENSEX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch
## $fitted.values): longer object length is not a multiple of shorter object
## length
rtns <- cbind(std.residual.DJIA[1067:4532],
              std.residual.GSPC[1067:4532],
              std.residual.NDX[2:3467],
              std.residual.GDAXI[953:4418],
              std.residual.FCHI[967:4432],
              std.residual.SSEC[658:4123],
              std.residual.SENSEX[674:4139]
              )

cor(rtns,method='pearson')
##             [,1]        [,2]          [,3]        [,4]          [,5]
## [1,]  1.00000000  0.97482107  0.4626361048 -0.02897670 -0.0448450319
## [2,]  0.97482107  1.00000000  0.4877918823 -0.02413000 -0.0454992144
## [3,]  0.46263610  0.48779188  1.0000000000 -0.05170464  0.0005284583
## [4,] -0.02897670 -0.02413000 -0.0517046389  1.00000000 -0.0142933609
## [5,] -0.04484503 -0.04549921  0.0005284583 -0.01429336  1.0000000000
## [6,] -0.02312062 -0.02512210  0.0025580560  0.01929367  0.0158593348
## [7,] -0.01925440 -0.02025158 -0.0048500626  0.01343810 -0.0123316220
##              [,6]         [,7]
## [1,] -0.023120619 -0.019254398
## [2,] -0.025122102 -0.020251582
## [3,]  0.002558056 -0.004850063
## [4,]  0.019293672  0.013438099
## [5,]  0.015859335 -0.012331622
## [6,]  1.000000000  0.013798025
## [7,]  0.013798025  1.000000000
cor(rtns,method='kendall')
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]  1.000000000  0.840936111  0.339267541 -0.018745030 -0.020497615
## [2,]  0.840936111  1.000000000  0.370060343 -0.012876103 -0.020819688
## [3,]  0.339267541  0.370060343  1.000000000 -0.026871801 -0.001039327
## [4,] -0.018745030 -0.012876103 -0.026871801  1.000000000 -0.012353858
## [5,] -0.020497615 -0.020819688 -0.001039327 -0.012353858  1.000000000
## [6,] -0.003618578 -0.004117842 -0.000861471  0.010674214 -0.002154094
## [7,] -0.004118841 -0.003936654 -0.004190783 -0.001609867 -0.020418012
##              [,6]         [,7]
## [1,] -0.003618578 -0.004118841
## [2,] -0.004117842 -0.003936654
## [3,] -0.000861471 -0.004190783
## [4,]  0.010674214 -0.001609867
## [5,] -0.002154094 -0.020418012
## [6,]  1.000000000  0.009381258
## [7,]  0.009381258  1.000000000
cor(rtns,method='spearman')
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]  1.000000000  0.962107436  0.454675399 -0.027960277 -0.030258828
## [2,]  0.962107436  1.000000000  0.484046514 -0.019486288 -0.030636594
## [3,]  0.454675399  0.484046514  1.000000000 -0.040064311 -0.001568855
## [4,] -0.027960277 -0.019486288 -0.040064311  1.000000000 -0.018482587
## [5,] -0.030258828 -0.030636594 -0.001568855 -0.018482587  1.000000000
## [6,] -0.005895352 -0.006491070 -0.001236738  0.015962569 -0.002794678
## [7,] -0.006206004 -0.006164209 -0.006189422 -0.002249416 -0.030567593
##              [,6]         [,7]
## [1,] -0.005895352 -0.006206004
## [2,] -0.006491070 -0.006164209
## [3,] -0.001236738 -0.006189422
## [4,]  0.015962569 -0.002249416
## [5,] -0.002794678 -0.030567593
## [6,]  1.000000000  0.013779081
## [7,]  0.013779081  1.000000000
#pairs(rtns)
splom2(rtns)

g <- rtns[,1]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    1   10   70  366 1183 1401  377   53    3    2
## 
## $density
##  [1] 0.0002885170 0.0028851702 0.0201961916 0.1055972302 0.3413156376
##  [6] 0.4042123485 0.1087709175 0.0152914022 0.0008655511 0.0005770340
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,2]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    4   14   82  362 1139 1409  398   52    4    2
## 
## $density
##  [1] 0.001154068 0.004039238 0.023658396 0.104443162 0.328620889
##  [6] 0.406520485 0.114829775 0.015002885 0.001154068 0.000577034
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,3]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
## 
## $counts
##  [1]    1    3   13  132  423 1072 1259  464   90    7    1    1
## 
## $density
##  [1] 0.0002885170 0.0008655511 0.0037507213 0.0380842470 0.1220427005
##  [6] 0.3092902481 0.3632429313 0.1338718984 0.0259665320 0.0020196192
## [11] 0.0002885170 0.0002885170
## 
## $mids
##  [1] -5.5 -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5  5.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,4]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    3   19  119  431 1051 1258  478   88   16    3
## 
## $density
##  [1] 0.0008655511 0.0054818234 0.0343335257 0.1243508367 0.3032313907
##  [6] 0.3629544143 0.1379111368 0.0253894980 0.0046162724 0.0008655511
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,5]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    2   20  122  422 1118 1205  479   82   10    6
## 
## $density
##  [1] 0.000577034 0.005770340 0.035199077 0.121754183 0.322562031
##  [6] 0.347663012 0.138199654 0.023658396 0.002885170 0.001731102
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,6]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6
## 
## $counts
##  [1]    1    1   10   42  112  424 1099 1166  460  123   22    4    2
## 
## $density
##  [1] 0.000288517 0.000288517 0.002885170 0.012117715 0.032313907
##  [6] 0.122331218 0.317080208 0.336410848 0.132717830 0.035487594
## [11] 0.006347374 0.001154068 0.000577034
## 
## $mids
##  [1] -6.5 -5.5 -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5  5.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
g <- rtns[,7]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall") 
    xfit<-seq(min(g),max(g),length=40) 
    yfit<-dnorm(xfit,mean=mean(g),sd=sd(g)) 
    yfit <- yfit*diff(h$mids[1:2])*length(g) 
    lines(xfit, yfit, col="black", lwd=2)

h
## $breaks
##  [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
## 
## $counts
##  [1]    3    3   65  373 1228 1346  396   48    3    1
## 
## $density
##  [1] 0.0008655511 0.0008655511 0.0187536065 0.1076168494 0.3542989036
##  [6] 0.3883439123 0.1142527409 0.0138488171 0.0008655511 0.0002885170
## 
## $mids
##  [1] -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5
## 
## $xname
## [1] "g"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

DJIA,GSPC

a <- rtns[,1]
b <- rtns[,2]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    7
## Name:  BB1
## 
## Parameter(s)
## ------------
## par:  1.05  (SE = 0.06)
## par2: 4.26  (SE = 0.1)
## 
## Dependence measures
## -------------------
## Kendall's tau:    0.85 (empirical = 0.84, p value < 0.01)
## Upper TD:         0.82 
## Lower TD:         0.86 
## 
## Fit statistics
## --------------
## logLik:  5119.14 
## AIC:    -10234.27 
## BIC:    -10221.97
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##     rho.1 
## 0.9714474
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##     rho.1        df 
## 0.9705497 4.1083364
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df = df),
                                      list(df = df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##   param 
## 7.36012
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2),margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

DJIA,NDX

a <- rtns[,1]
b <- rtns[,3]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  0.54  (SE = 0.01)
## par2: 2  (SE = 0.11)
## 
## Dependence measures
## -------------------
## Kendall's tau:    0.37 (empirical = 0.34, p value < 0.01)
## Upper TD:         0.42 
## Lower TD:         0.42 
## 
## Fit statistics
## --------------
## logLik:  808.31 
## AIC:    -1612.62 
## BIC:    -1600.31
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##     rho.1 
## 0.4660699
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##     rho.1        df 
## 0.5363481 1.7946408
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df=df),
                                      list(df=df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##    param 
## 0.836929
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))



sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

DJIA,GDAXI

a <- rtns[,1]
b <- rtns[,4]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  -0.03  (SE = 0.02)
## par2: 13.28  (SE = 3.49)
## 
## Dependence measures
## -------------------
## Kendall's tau:    -0.02 (empirical = -0.02, p value = 0.1)
## Upper TD:         0 
## Lower TD:         0 
## 
## Fit statistics
## --------------
## logLik:  9.75 
## AIC:    -15.5 
## BIC:    -3.19
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##       rho.1 
## -0.03054467
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##      rho.1         df 
## -0.0300592 13.2810415
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df=df),
                                      list(df=df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##       param 
## -0.01112868
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))



sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

DJIA,FCHI

a <- rtns[,1]
b <- rtns[,5]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  -0.03  (SE = 0.02)
## par2: 9.21  (SE = 1.66)
## 
## Dependence measures
## -------------------
## Kendall's tau:    -0.02 (empirical = -0.02, p value = 0.07)
## Upper TD:         0.01 
## Lower TD:         0.01 
## 
## Fit statistics
## --------------
## logLik:  21.35 
## AIC:    -38.7 
## BIC:    -26.4
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##       rho.1 
## -0.03932221
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##       rho.1          df 
## -0.03382907  9.20959786
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df=df),
                                      list(df=df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##       param 
## -0.01643638
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))



sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

DJIA,SSEC

a <- rtns[,1]
b <- rtns[,6]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  -0.01  (SE = 0.02)
## par2: 16.17  (SE = 5.01)
## 
## Dependence measures
## -------------------
## Kendall's tau:    -0.01 (empirical = 0, p value = 0.75)
## Upper TD:         0 
## Lower TD:         0 
## 
## Fit statistics
## --------------
## logLik:  6.34 
## AIC:    -8.67 
## BIC:    3.63
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##       rho.1 
## -0.01696391
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##      rho.1         df 
## -0.0129055 16.1649914
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df=df),
                                      list(df=df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##       param 
## -0.01336535
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))



sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

DJIA,SENSEX

a <- rtns[,1]
b <- rtns[,7]

m <- pobs(as.matrix(cbind(a,b)))

u <- m[,1]
v <- m[,2]

mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)

plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)

selectedCopula <- BiCopSelect(u,v,familyset=NA)

summary(selectedCopula)
## Family
## ------ 
## No:    124
## Name:  Rotated Tawn type 1 90 degrees
## 
## Parameter(s)
## ------------
## par:  -1.1  (SE = 0.08)
## par2: 0.06  (SE = 0.05)
## 
## Dependence measures
## -------------------
## Kendall's tau:    -0.01 (empirical = 0, p value = 0.72)
## Upper TD:         0 
## Lower TD:         0 
## 
## Fit statistics
## --------------
## logLik:  4.14 
## AIC:    -4.28 
## BIC:    8.02
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
##       rho.1 
## -0.01430268
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)

fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
##       rho.1          df 
## -0.01250382 34.12112734
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
                    paramMargins=list(list(df=df),
                                      list(df=df)))

sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)

fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
##       param 
## -0.01099794
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))



sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
##    param 
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mu.1, sd=sd.1),
                                      list(mean=mu.2, sd=sd.2)))

sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)